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Resumen: Se emplea una ecuación de cable no lineal homogeneizada para describir la 
interacción entre un campo eléctrico externo y una fibra nerviosa mielínica. Aplicando el 
concepto de intervalo de influencia del campo sobre el axón, y un modelo simplificado de 
membrana excitable, se utiliza un método de Galerkin para aproximar la ecuación a deri- 
vadas parciales no lineal del modelo original por un sistema de ecuaciones diferenciales 
ordinarias en las amplitudes de modo resultantes. Se presentan algunos resultados de 
simulación digital en un espacio de amplitudes de modo, para casos en los que la variable 
de recuperación se fija en su valor de reposo. Se identifican una amplitud dominante y se 
utiliza para simplificar el sistema de ecuaciones diferenciales ordinarias. Se estudia la 
dinámica del umbral en el sistema de orden reducido y se construye un modelo matemáti- 
co de las curvas intensidad-duración aplicable a la estimulación eléctrica y magnética. Se 
obtienen fórmulas analíticas para la curva intensidad-duración, en particular para la cro- 
naxia y la reobase correspondientes a la estimulación catódica de cierre. Se analizan la 
estimulación anódica de cierre y se sugiere como abordar un análisis del bloqueo anódico 
por el cátodo. Se ponderan las limitaciones del modelo de axón y de las técnicas de análi- 
sis utilizadas. 

Palabras clave: Ecuación del cable homogeneizada, curvas intensidad-duración, umbra- 
les, cronaxia, reobase, estimulación eléctrica funcional, estimulación magnética, fibra 
nerviosa mielínica, análisis modal no lineal, bifurcaciones globales, estimulación catódica 
de cierre, estimulación anódica de cierre. 


Abstract: A nonlinear and homogenized cable equation is employed to describe the in- 
teraction between an external electric field and a myelin nervous fiber. Applying the con- 
cept of the interval of influence of the field over the axon, and a simplified model of ex- 
citable membrane, Galerkin?s method is used to approximate the nonlinear partial differ- 
ential equation of the original model by a system of nonlinear differential equations in the 
resulting mode amplitudes. Some digital simulation results in the space of mode ampli- 
tudes are reported, for the recovery variable fixed at its rest value. Dominant mode ampli- 
tude is identified and used to simplify the system of ordinary differential equations. 
Threshold dynamics is studied in the reduced order system and a mathematical model of 
strength-duration curves for electric and magnetic stimulation is constructed. Analytical 
formulae for strength-duration curves, in particular for chronaxie and rheobase, are ob- 
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tained for cathodic make excitation. Anodic make excitation from the cathode is dis- 
cussed. An approach to anodic block from the cathode is suggested. The limitations of the 
axon model and the analytical techniques employed are assessed. 

Key words: Homogenized cable equation, strength-duration curves, thresholds, chronax- 
ie, rheobase, functional electric stimulation, magnetic stimulation, myelin nervous fiber, 
non-linear modal analysis, global bifurcations, cathodic make excitation, anodic make 
excitation. 


(1)Introducción 


El objetivo de la estimulación eléctrica y magnética de células excitables con 
fines terapeúticos, es generar o bloquear potenciales de acción en sitios adecua- 
dos a los fines perseguidos [1], [14], [20]. En ambos casos el mecanismo de esti- 
mulación involucra una interacción entre un campo eléctrico externo y la célula 
excitable que se desea estimular. 

En el caso de la estimulación eléctrica el campo se produce mediante un arreglo 
de electrodos, ubicados más o menos próximos al tejido cuya actividad se preten- 
de estimular o bloquear. A través de los electrodos se inyecta una corriente eléc- 
trica en el conductor de volumen formado por los tejidos biológicos, con una 
duración y una amplitud que se seleccionan buscando que la estimulación resulte 
exitosa y específica. 

En el caso de la estimulación magnética, la variación temporal del campo magné- 
tico genera un campo eléctrico inducido, y es este último el responsable de la 
estimulación o el bloqueo en la célula blanco. El flujo magnético variable se logra 
inyectando una corriente eléctrica variable en una bobina o conjunto de bobinas. 
La variación del flujo magnético en el conductor de volumen genera una distribu- 
ción espacial de corrientes y campos eléctricos inducidos. 

El proceso de estimulación eléctrica o magnética en el sistema formado por el 
dispositivo estimulador, el conductor de volumen de los tejidos y el elemento o 
elementos que se desea estimular, puede describirse a nivel de caja negra como 
una correspondencia entre cada historia de corriente eléctrica aplicada (inyectada 
en el conductor de volumen de los tejidos mediante electrodos o inyectada en las 
bobinas para producir el flujo magnético variable) y una variable binaria A que 
toma el valor 0 si la estimulación fracasa y 1 si tiene éxito [15], [19]. 

La variable binaria refleja el carácter de fenómeno de umbral que presenta la 
estimulación eléctrica y magnética. Se puede obtener a partir de medidas eléctri- 
cas u observando manifestaciones externas como contracciones musculares o 
inhibición de funciones. El esquema de la Fig. 1 resume el umbral de estimula- 
ción eléctrica o magnética desde el punto de vista de una caja negra. 
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Fig. 1. La estimulación eléctrica y magnética considerada como caja negra. 


La determinación de la curva intensidad-duración es una herramienta que se utili- 
za en la práctica para llevar a cabo la identificación del umbral definido mediante 
la variable A. 

Utilizando pulsos rectangulares de corriente controlada, que generan un pulso de 
campo eléctrico aproximadamente rectangular en los tejidos biológicos, hay dos 
formas de determinar la curva de intensidad-duración en el caso de la estimula- 
ción eléctrica: 

a) se fija la duración del pulso estimulador y se varía su amplitud. 

b) se fija la amplitud y se varía la duración. 

En ambos casos se determina el par amplitud-duración umbral'. Los resultados 
obtenidos se pueden graficar en un plano, ubicando la duración en las abscisas y 
la intensidad en las ordenadas. Los pares duración-amplitud umbral forman una 
curva que se asemeja a una hipérbola y que separa en primer cuadrante en dos 
regiones: una región de estímulos sub-umbrales, debajo de la curva, y una región 


| Se tienen en cuenta pequeñas variaciones aleatorias introduciendo un criterio estadístico 
para determinar el umbral. Por ejemplo, se considera pulso umbral cuando la excitación 
resulta exitosa en el 50 % de los casos para una misma amplitud y duración pre-fijadas en 
el equipo estimulador [9]. 
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de estímulos supra-umbrales, encima de la curva. 

Con los equipos para estimulación magnética de tercera generación se pueden 
inyectar pulsos de corriente en las bobinas que generan en los tejidos pulsos casi 
rectangulares de campo eléctrico inducido [11]. Esto permite determinar curvas 
intensidad-duración en forma análoga a la empleada en la estimulación eléctrica. 
La determinación del umbral de estimulación eléctrica o magnética debe hacerse 
en cada caso, porque la curva intensidad-duración depende de la posición y tama- 
ño de los electrodos o las bobinas, de las propiedades del conductor de volumen 
formado por los tejidos biológicos y de la naturaleza y posición de la o las células 
blanco de la estimulación. 


El propósito principal del presente trabajo es construir un modelo matemático 
simplificado de las curvas intensidad duración cuando la estimulación se produce 
en una fibra nerviosa mielínica. 

La Fig. 2 muestra un esquema de una porción de un axón envuelto en manguitos 
de mielina aportadas por células de Schwamn. A intervalos regulares la envoltura 
aislante de mielina se interrumpe y la membrana del axón queda al descubierto, 
formando los nodos de Ranvier, cuyas membranas son excitables. 
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Fig. 2 Esquema de un tramo de fibra nerviosa mielínica. 


Se indica un electrodo puntual ubicado por encima y a cierta distancia de un no- 
do, a través del cual se inyecta una corriente / . Con v 


e,n—1 


Von y Ve n+1 se repre- 


senta los potenciales eléctricos en el exterior de tres nodos consecutivos. La lon- 
gitud entre nodos es £, la longitud de un nodo es y, su diámetro es d y el diá- 


metro externo del manguito de mielina es D . 


Se representa el axón mediante lo que se suele denominar una ecuación del cable 
homogeneizada, como la obtenida por Basser con el propósito de describir el 
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comportamiento de la fibra por debajo del umbral [1]. Aquí se extiende su apli- 
cación para describir perturbaciones de las variables de estado de la membrana 
hasta alcanzar un estado umbral. La membrana a la que se hace referencia es una 
membrana equivalente, cuyas variables de estado aparecen como variables de- 
pendientes del tiempo y la posición en la ecuación homogeneizada. Estas varia- 
bles de estado emergentes pueden interpretarse como una especie de promedio de 
las correspondientes variables en los nodos de Ranvier y en los manguitos de 
mielina. 

Se formula la ecuación del cable homogeneizada para un modelo de membrana 
excitable del tipo de Fitz Hugh-Nagumo, que emplea dos variables de estado: el 
voltaje transmembrana y una variable de recuperación [7]. 

En esta primera parte del trabajo se estudia la dinámica del umbral para la varia- 
ble de recuperación fija en su valor de reposo. En la segunda parte se estudiará la 
dinámica del umbral cuando la variable de recuperación se encuentra activa y se 
compararán los resultados con los obtenidos cuando se encuentra fija. 


(2)Modelo matemático de la interacción campo-axón basado en una 
ecuación de cable no lineal homogeneizado 


La figura 3 muestra un tramo de un axón mielínico, curvado e inmerso en un 
campo eléctrico externo. En este caso se supuso que el campo externo es unifor- 
me, pero en general puede variar en magnitud y dirección de forma regular pero 
por lo demás arbitraria. 


? 


sentido 
positivo 


Fig.3 Tramo de axón curvado en un campo eléctrico uniforme. Se advierte que la com- 
ponente del campo tangente al axón se modifica de un nodo de Ranvier a otro. 
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(2.1) Sistema de ecuaciones a derivadas parciales no lineal para el cable ho- 
mogeneizado 

Un modelo detallado requiere emplear una ecuación de cable separada para cada 
nodo y otra para cada manguito de mielina, introduciendo condiciones en cada 
interface nodo-manguito para compatibilizar las soluciones entre sí [7], [8]. 

No obstante, mediante un proceso de homogeneización, basado en el método 
asintótico de las escalas múltiples (dos escalas espaciales en este caso, una del 
orden de la longitud del nodo y otra del orden de la distancia entre nodos) se pue- 
de hallar una ecuación de cable emergente, como una primera aproximación, a 
partir de operar sobre un modelo de la micro-estructura detallada del axón [1]. 
Para los fines del presente trabajo se puede partir del siguiente sistema: 


A DE 
A a ea e, par (la) 
Ot L OS OS 
ow pu 
T,—=—lvV-Yy:w 1b 
"Ot E 70) (1b) 


En las ecuaciones (la) y (1b) el campo escalar vlt, s) representa el potencial a 
través de la membrana excitable equivalente (es decir, emergente del proceso de 
homogeneizado), relativo a su valor de reposo, en el instante £ y en la posición 
de abscisa curvilínea sa lo largo del axón (se introduce una abscisa curvilínea 
para poder trabajar con fibras curvadas). El campo escalar wlt, s) representa una 
variable de recuperación de la membrana equivalente. 
El término que representa las no linealidades que hacen posible la excitación y el 
efecto de la variable de recuperación viene dado por: 

Nlv, w)=bv? —cv' —a-w (2) 
Corresponde a un modelo simplificado de membrana excitable del tipo conocido 
como de Fitz-Hug y Nagumo [7]. Los parámetros que aparecen en estas ecuacio- 
nes son positivos. 
Por definición: > = ( 2) . + e e (€) 

A LA LA 


mie 
A es una constante de espacio del cable homogeneizado, expresada en función de 


las constantes de espacio 4,,,, del manguito de mielina y 4,,, del nodo de Ranvier. 


uN RE lA 
También por definición: | A (4) 
L Amic L Ano 
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T es una constante de tiempo del cable homogeneizado, expresada en función de 
las constantes de tiempo 7,,, del manguito de mielina y 7,,, del nodo de Ranvier. 
El término que conecta el estado de la membrana equivalente con el campo eléc- 
trico externo, involucra la derivada de la componente del campo externo tangen- 
cial al axón respecto de la abscisa curvilínea del axón. En las condiciones usuales 
de estimulación eléctrica y magnética se puede expresar como un producto de una 
función del tiempo por una función del espacio (de la abscisa curvilínea): 


- Pato =F(s): g(t) (5) 


En la estimulación eléctrica gls) es la corriente inyectada por el electrodo activo 


i el), mientras que en la estimulación magnética es la derivada respecto del 


tiempo de la corriente en la bobina = la (e). 
Í 


En ambos casos F (s) representa la variación espacial de la componente del 
campo externo paralela al axón, por unidad de longitud tomada a lo largo de la 
fibra. Lo denominaremos factor de forma. A vía de ejemplo, para una fibra 
recta y un electrodo puntual situado a una distancia z de la fibra, en un 
punto de abscisa (tomada sobre la fibra) x=0, el factor de forma, cuando 
el conductor de volumen es homogéneo e isótropo y posee una conductivi- 
dad «o, viene dado por la fórmula: 


1 2x*-z 


F(x)= ano (Y (6) 


X 


(2.2) Hipótesis del intervalo de influencia del campo externo sobre el cable 
homogeneizado y deducción de las ecuaciones para las amplitudes de modo. 
En lo que sigue introduciremos la hipótesis del intervalo de influencia? [15] [16]: 


? Esta suposición es análoga a la hipótesis de las longitudes extrapoladas que se efectúa en 
la física de reactores nucleares [17]. Permite describir los campos que caracterizan el 
estado de la membrana como suma de términos construidos a partir de las funciones pro- 
pias (modos) de un operador diferencial lineal espacial (operador de Laplace), en forma 
análoga a como se expresan el flujo neutrónico y otros campos en el núcleo de un reactor. 
Esos últimos se suponen nulos sobre la frontera extrapolada del reactor, lo cual no es ni 
más ni menos que una buena aproximación (puesto que siempre se puede detectar algún 
neutrón, proveniente del reactor, más allá de la frontera extrapolada). 
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admitiremos que (s) se anula fuera de un intervalo acotado de longitud / . 

Esta suposición se basa en la constatación de que, por lo general, el proceso de 
estimulación, desde el reposo y hasta que se alcanza el umbral es un fenómeno 
espacialmente localizado: fuera de una región de influencia acotada es posible 
despreciar la perturbación que el campo externo produce en las variables de esta- 
do de la membrana y asumir que estas últimas permanecen fijas en sus valores de 
reposo. En el caso de los campos vít,s) y wlt,s) se supone que se anulan en los 


extremos del intervalo de influencia. Si el origen se pone en el punto medio de 


dicho intervalo, se supone entonces que: da 5) =0 y 91 3 z=0 


que aparece en la ecuación 


01...) 
OS 


Consideremos ahora el operador £|...]= —E z 


(la) actuando sobre vlt,s). Planteamos el problema de hallar las funciones pls) 


2 
tales que -lols)l- o-p(s) siendo pls) nula en los extremos del intervalo 
s 


de influencia de los electrodos o las bobinas sobre la fibra. Con el origen de coor- 
denadas en el punto medio del intervalo de influencia, se obtienen los siguientes 


Zi Z 
valores propios: (0, = — siendo £ la longitud del intervalo de influencia y 
n=1, 2, 3,... Cuando n es impar 0)= 5025), y cuando es par 
Ñ Ed L 


p, (s)= (ise 122) Estas funciones propias están normalizadas de modo que 


pleto que permite desarrollar en serie de Fourier cualquier función regular en el 

intervalo - A E con suficiente aproximación si se emplea un número de tér- 
Zi *2 

minos lo bastante grande. Teniendo esto último en cuenta planteamos desarrollos 


aproximados en términos de las funciones propias dependientes de la posición 
sobre el cable, con amplitudes de modo dependientes del tiempo: 
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5(8s)=Y 4,(0-0,(s) (13) Ea) y Bl0:0(s) 05) 


Ahora reescribimos las ecuaciones (la) y (1b) bajo la forma 


Ov 17, Ov 
E VE SAA -Els)- g(t)=R, [v, w]=0 (8a) 
10) 
Ea 7-1 w)=R,[v,w]=0 (8b) 


Ahora aplicamos el método de Galerkin [4], sustituyendo (7a) y (7b) en (8a) y 
(8b) y proyectando sobre las primeras M funciones propias e igualando a cero 


esas proyecciones: 
Ll 


fo, (s)-R, [9,7]. ds =0 (9a) p,(s)-R,[9,%]-ds =0 (9b) 


OS 


Se obtiene un sistema de ecuaciones diferenciales ordinarias no lineales para las 
amplitudes de modo: 


dA, A un Y Y 
a mn? YA? 1022 ap Ayo 2 Sapos A A, A 0:B,+4 E, g(t) 
(10a) 
E o =E(a, 85) (n=1, 2, 3,.... M) (10b) 


1 
2 
En estas expresiones se define  0,,,= j p,[s)-p,(s)- p,(s)- ds , Mientras que 


npgr 


a 
ps j p,[(s)-9,[s)-p,(s)- p,(s)- dx Estos coeficientes dependen solamente de las 
l 


Sl] 


funciones propias del operador: son independientes de la naturaleza del sistema 
formado por el dispositivo estimulador, el conductor de volumen y la fibra blanco 
de la excitación. Sus expresiones analíticas pueden obtenerse mediante el uso de 
paquetes de software para cálculo simbólico o recurriendo a tablas de integrales. 
En la referencia [10] aparecen los resultados para los índices de modo com- 
prendidos entre n=1 y n=09. 
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Las proyecciones de los factores de forma vienen dados por las fórmulas? 


F, = | p,(s)-F(s)- ds (n=1, 2, 3,.... M) (11) 


El sistema de ecuaciones que incluye una ecuación a derivadas parciales no li- 
neales y describe la dinámica de la membrana del cable homogeneizado queda así 
sustituido por un sistema de ecuaciones diferenciales ordinarias para las amplitu- 
des de modo 4, (e) y B,(t) correspondientes a las expresiones aproximadas 
(fórmulas (7a) y (7b)) de las variables de estado. Cada ecuación de evolución de 
un modo presenta un término forzante proporcional a la corriente en el electrodo 
en el caso de la estimulación eléctrica, o a la velocidad de variación de la co- 
rriente en la bobina en el caso de la estimulación magnética. 

En el factor de proporcionalidad aparece la proyección del factor de forma dada 
por la ecuación (11). Si esta proyección es nula o al menos despreciable, el 
modo solo puede activarse a través de su acoplamiento con los otros mo- 
dos, algunos de los cuales se verán forzados por unas proyecciones FF, no 


nulas. 
Aplicando el lema de Riemann-Lebesgue [3] a la ecuación (11) se deduce 
que cuando el índice de modo n tiende a infinito, FF, tiende a cero. 


La Fig. 4 muestra un intervalo de influencia. 


21 FO) 


400 | 500 800 400 500 600 400 500 $00 
Intervalo de influencia 
Fig. 4 El factor de forma para la estimulación mediante un electrodo puntual, el intervalo 
de influencia correspondiente y la aproximación en serie de Fourier utilizando los cinco 
primeros modos. Las componentes correspondientes a los modos de índice par se anulan. 


3 Cuando la estimulación se debe a un electrodo puntual y la fibra no se encuentra curva- 


da, las proyecciones del factor de forma resultan proporcionales a 1/ E, 


Rev. SCP, Vol. 17, N” 2, pp. 129-154, Diciembre 2012 ISSN 0379-9123 


El sector izquierdo de la fig. 4 muestra una gráfica el factor de forma para 
un electrodo estimulador puntual, dado por la ecuación (6) (con a =1), y 
la construcción del intervalo de influencia para este factor de forma. El 
método de construcción de los intervalos de influencia para electrodos de 
formas y tamaños cualesquiera y para arreglos de electrodos puede hallarse 
en las referencias [15] y [17]. 

En el sector central de la figura se representan las primeras tres componen- 
tes no nulas del desarrollo de Fourier del factor de forma considerado: 


F, cos dea ES y ES 
L L E 


En el sector derecho de la figura se puede ver la reconstrucción del factor 
de forma como suma de esas tres componentes, de índice impar (pero que 
en cuanto funciones de x, son funciones pares): 


Ecos 2 + F,cos| == 2 +F¿0 E re 
E 7 L 


(3)Análisis modal no lineal de la ecuación del cable homogeneizada: 
dinámica del umbral para una membrana con la variable de recupe- 
ración fija. 


Supongamos ahora que la variable de recuperación permanece fija en su 
valor de reposo. e ese caso las ecuaciones (10) se reducen al sistema: 
dA, 


a E a +£o y Bo. Ay A, ek y O04 Ay Ay A+ AF, gt) 


pq=1 pqur=1 


(12) 
(3.1) Simulación digital de la estimulación catódica de cierre. 
La Fig. 5 muestra los resultados de una simulación digital? efectuada a 
partir de las ecuaciones (12) para M =5, con el factor de forma simétrico 
del electrodo puntual, empleando la descomposición que muestra la Fig. 4. 
En este caso el campo eléctrico externo activa solo tres amplitudes de mo- 
do (4,, 4, y 45). La dinámica de la fibra se describe entonces en un 


espacio de estados de tres dimensiones. Para hacer la simulación digital se 


1 Las corridas de simulación digital fueron realizadas por el ingeniero Héctor Korenko, en 
Matlab y utilizando los valores de los parámetros indicados por el autor. 
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tomó: 4p=1 4e=1 7+=3ms A4=1mm y una distancia electrodo fibra de 
L L 8 


z=20mm. 


Amplitud de modo 


% 0.5 1 


tiempo (ms) tiempo (ms) 
Fig. 5: Dinámica forzada por un pulso catódico y dinámica libre con tres amplitudes de 
modo. 


Para z=20mm se obtiene: F¡ =-9,5x107 F,=-2,4x10* y F,=-1,6x10% 
Como gls) se eligió un pulso rectangular de amplitud negativa, que en 
este modelo corresponde a un pulso catódico. La duración del pulso se 
fijó en 1 ms, de modo que resultara corto respecto de la constante de tiem- 
po de la membrana de la fibra homogeneizada. 

La figura de la izquierda muestra la activación de las amplitudes de los tres 
modos de índice impar a partir de sus valores de reposo, mientras dura el 
pulso estimulador. 

La figura de la derecha muestra la evolución de esas mismas amplitudes de 
modo una vez que el pulso estimulador ha finalizado. 

Se ve que el modo fundamental, de amplitud 4,, es dominante respecto 
de los otros dos modos, pese a que al final del pulso estimulador su ampli- 
tud es inferior a las amplitudes de los otros dos modos. 

En el caso que se muestra en la Fig. 5, la amplitud del pulso es lo bastante 
grande como para que, una vez finalizado el pulso estimulador, 4 , tienda 


asintóticamente a un valor no nulo. Durante su evolución, 4, arrastra a las 


otras amplitudes, que se mantienen en cuasi-equilibrio con la amplitud de 
modo dominante. 
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Los valores asintóticos de las amplitudes de modo corresponden a un esta- 
do excitado de la fibra, espacialmente no uniforme y estable (debido a que 
la variable de recuperación se mantiene en su valor de reposo). 

El pulso estimulador resulta supra-umbral. Si se disminuye lo suficiente la 
amplitud del pulso, una vez finalizado éste las tres amplitudes de modo 
regresan a sus valores de reposo: el pulso estimulador resulta sub-umbral. 
Si se alarga la duración del pulso estimulador por encima de 2 ms, al final 
del pulso 4 resulta mayor que 4, y esta última amplitud resulta mayor 


que A5. Pero luego de finalizado el pulso el primer modo domina y arras- 


tra a los demás, igual que ocurre para pulsos cortos, ya sea hacia el reposo 
o hacia el estado excitado. 


(3.2) Dinámica del umbral para el modo dominante desacoplado. 

Estos y otros resultados análogos de simulación digital sugieren que el 
primer modo resulta ser dominante cuando la estimulación se produce me- 
diante pulsos catódicos. 

Admitiendo que esto es así, en una primera aproximación se lo puede des- 
acoplar de las otras amplitudes de modo. 

Con esta simplificación obtenemos una dinámica en un espacio de estados 
unidimensional que en principio podría aplicarse tanto a la estimulación 
eléctrica como magnética, si el intervalo de influencia y la proyección del 
factor de forma se interpretan adecuadamente: 


dA, 


q 
a =-(1+7? A +E OO Aj CO A +1? E g(0) (13) 


Los parámetros O valen: 0,,, =1,2 y Q,,,, =1,5 


Una vez finalizado el pulso estimulador, 4, queda abandonada a sí mis- 
ma, con una dinámica determinada por la ecuación: 


A 2 
CA de A+ L0oÓ Af 008 A (14) 
Los puntos de equilibrio posibles son las raíces reales de 
p q p 
A 2 
-(1+11? ad +2b-8, A/ 08 Aj =0 


Una raíz real, para todos los valores de los parámetros, es 4,=0, que co- 


rresponde al estado de reposo de la membrana de la fibra. 
Las otras dos raíces, posiblemente complejas, son: 
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b-0 b-0 
Lu m JA 1 a) A, = a, +/A (Ss b) 
208 111 2001111 


El discriminante A viene a su vez dado por la fórmula: 
2 2 
[29m O (17) 
20 01111 E 0 ) e 
LL 1111 
Cuando A es positivo, además del estado de reposo, hay otros dos estados 


de equilibrio dados por las fórmulas (15 a), (15 b) y (16). 
En este caso A, , es un estado de equilibrio instable: si la amplitud inicial 


se encuentra entre 4,=0 y A,,, la amplitud del modo dominante retorna 


siempre al estado de reposo. A este intervalo lo denominaremos conjunto 
de decaimiento. 
Si la amplitud inicial se encuentra comprendida entre A, ,, y 4,,la ampli- 


Lu 
tud del modo dominante tiende siempre hacia el estado excitado 4,.,. A 
este intervalo lo denominaremos conjunto amplificador. 

Por tanto, A,, es una amplitud umbral que debe ser superada para que la 


excitación de la fibra pueda producirse. 
La Fig. 6 muestra un diagrama de bifurcación global para los puntos de 
equilibrio de la amplitud del modo dominante. 


dan ] 


Lu 


1.00 — 


17) 


l 
Fig. 6: Diagrama de bifurcaciones en función de un parámetro P G 0—) 
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20 Al Edic 
Cuando el parámetro de bifurcación P (es decir a O E) es inferior a un 


valor crítico P,,, existe un único punto de equilibrio, el estado de reposo: 


cr > 
la fibra no es excitable. Pero si el parámetro de bifurcación es superior a su 
valor crítico, aparecen un umbral de excitación y un estado excitado. 


Para el valor crítico del parámetro de bifurcación se verifica A=0. 


; : l 
De la ecuación (16) se desprende que los valores críticos de ” y de E se 


> b.a? 0 
relacionan entre sí a través de ll (£ ) Sis ”/ () (18) 
4c08 11 L A 


(4) Deducción de una fórmula analítica para la curva intensidad dura- 
ción en la excitación catódica de cierre con un modo dominante. 


Linealizando la ecuación (13) en el conjunto de decaimiento se obtiene: 


2 qe 
o ES Jaca, (19) 


Ahora, consideramos un pulso g(t) que comienza en t=0 con la mem- 


brana en reposo: 4,=0 
Supongamos que al finalizar el pulso, de duración d , se verifica 
A(Y)=A,, (20) 


Entonces ese pulso es un pulso umbral. 
Integrando la ecuación (19) en estas condiciones se obtiene, al finalizar el 


EN? t-t 

pulso: A (d)= (25) J 0 el . e(t). dt' Q1) 
T 0 Ss 

De (0) y (1) resulta que el pulso será un pulso umbral siempre y cuando 


0 
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En esta fórmula t, es la constante de tiempo del sistema formado por el 


dispositivo estimulador, el conductor de volumen y la fibra excitable, y 


viene dada por fa (Q3) 


S 77? de: 2 
1+ 7 
Restrinjamos la estimulación considerada a la estimulación eléctrica me- 
diante un cátodo activo. Si el pulso es un pulso rectangular de corriente 
catódica de amplitud / , de la condición (22) se deduce que el pulso cató- 
dico será umbral si verifica la siguiente relación entre su amplitud y su 
duración: 


js (24) 


= ep - a 


En esta fórmula R es la corriente reobásica o reobase, por debajo de la 
cual la fibra no puede ser excitada, no importa cuán larga sea la duración 


a A 


7 lu (25) 


del pulso estimulador: R=|1+ 5 
E 
La cronaxia t,, es la duración del pulso para la cual la amplitud umbral es 


el doble de la reobase. De acuerdo con la fórmula (24) es proporcional a la 
constante de tiempo del sistema: t., 0.63 -£. (26) 


(5)Simulación digital de la excitación anódica de cierre. 


Las ecuaciones de evolución de los modos de índices 3 y 5, desacoplados 
de los demás modos pero forzados por el dispositivo estimulador, vienen 
dadas por: 


qe 
a =-(1+9-7? DA + EDO A? ESO A? +2 E g(0) (Q7) 


q? 
A As HE bss AS Lor AS +4 Kg (8) 
En estas ecuaciones F, y Fj son, como siempre, negativos, 


33, = 055; =1,5 pero hay un cambio de signo: ahora 6,,, =—0,4 mientras 
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que 6... =+0,2 es positivo al igual que el valor ya visto para el modo uno: 
0, =>+1,2 

Debido al signo negativo de 0,,,, el modo tres desacoplado presenta un 
umbral 4, que (cuando existe) es negativo. El signo del umbral 4,,, 


para el modo 5 (cuando este umbral existe) es positivo al igual que el signo 
del umbral del primer modo. 

Ahora bien, si la amplitud de un pulso estimulador es positiva, es decir, si 
es un pulso anódico, y la membrana parte del estado de reposo, cabe espe- 
rar que las amplitudes de modo evolucionen alejándose del umbral 4, del 


modo uno desacoplado, aproximándose al umbral 4;,, del modo 3. 


Los resultados de la simulación digital con los tres modos acoplados con- 
firman esa presunción y muestran en el plano de las amplitudes de modo 
A,- Ay una curva umbral que puede ser traspasada tanto por órbitas aso- 


ciadas a pulsos catódicos como anódicos. Ver Fig. 7 


15 a T T 7 JE 
z = 5mm 
P X= ¿o 7 
1 b LC tc LC 1/50 
Electrodo puntual 
0.5f 1 
2 ob 7 
-0.5f 
Separatriz umbral 
AF | 
| 
1 L 
145 A 0,5 0 0.5 1 15 


A1 
Fig. 7. Separatriz umbral y órbitas asociadas a pulsos catódicos y anódicos. 


Se pueden observar dos haces curvos de órbitas que parten del origen, cal- 
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culados para pulsos de diferente duración e intensidad. 

El haz que asciende hacia la derecha a partir del origen corresponde a la 
estimulación catódica. Si una de esas curvas atraviesa la curva umbral, se 
produce una excitación catódica de cierre (E.C.C.). Si no la atraviesa, el 
estado vuelve al origen y la excitación habrá fracasado. 

Las curvas del haz que desciende hacia la derecha a partir de origen co- 
rresponden a la estimulación anódica. Si traspasan la curva umbral se pro- 
duce una excitación anódica de cierre (E.A.C.). Si no la traspasa, el esta- 
do vuelve a su condición de reposo. 

La curva umbral de la Fig. 7 es la traza en el plano 4,-4,de una superfi- 


cie umbral en el espacio 4,-4,-4As. 


(6)Discusión y conclusiones 


(a)El cociente entre la longitud del nodo de Ranvier de la fibra mielinizada 
y la longitud entre nodos verifica O < = <1. Cuando E =1se tiene el caso 
extremo de una fibra completamente desprovista de envoltura de mielina. 
De modo que poniendo 7 =1 en las fórmulas deducidas en este trabajo se 


obtienen las fórmulas correspondientes a las fibras amielínicas. Esto es una 
consecuencia de haber trabajado a partir de una ecuación de cable homo- 


geneizada. En el otro extremo, cuando 7= 0, los nodos han desaparecido 
y la fibra, completamente recubierta de mielina, ya no resulta excitable. 

De la fórmula (15 a) y de la Fig.6 se desprende que si 7 disminuye, el 
umbral para el modo uno desacoplado aumenta hasta un punto en el que 
coincide con el estado excitado y si 07 continúa disminuyendo, la fibra se 


vuelve no excitable. A partir de la relación (18) se puede hallar una expre- 


4c0, LV 
sión para ese valor crítico: (2 ) Sie o de (0) 
dos 0 Gr A 
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(b)La ecuación (19), obtenida por linealización, en el conjunto de decai- 
miento, de la ecuación de la dinámica del umbral para el modo uno des- 
acoplado, no es otra cosa que una versión del modelo de un factor, debido 
a Blair [7] [15] [17], para la estimulación catódica de cierre. 

Pero el resultado es muy diferente de ese modelo clásico. 

Por un lado, el modelo de Blair aparece extendido a la estimulación mag- 
nética. Por otro lado, en lugar de conceptos no bien definidos y variables 
fenomenológicas introducidas ad-hoc (“variable de excitación”, “umbral” 
y “factor de tiempo para la excitación”) que aparecen en el modelo de 
Blair original, en nuestro caso la variable de estado y los parámetros que 
aparecen poseen un significado electrofisiológico bien definido y se pue- 
den relacionar mediante fórmulas analíticas con las propiedades del siste- 
ma formado por el dispositivo estimulador, el conductor de volumen for- 
mado por los tejidos y la fibra cuya actividad eléctrica se busca excitar o 
bloquear. 


(c)Siempre que el electrodo o la bobina no se encuentren muy cerca de la 
fibra mielínica, cabe esperar que las fórmulas para la reobase (25) y la 
constante de tiempo (23) resulten aplicables para estudiar el efecto de la 
variación de los parámetros del sistema dispositivo estimulador-conductor 
de volumen-fibra blanco de la estimulación. La constante de tiempo de- 


pende de 7 través de T y de A. Cuando 7 aumenta de 0 a 1, la cons- 


tante de tiempo varía desde la correspondiente a una envoltura completa de 
mielina hasta la de una fibra completamente desprovista de mielina. 

Los incrementos en la constante de tiempo con la longitud del intervalo de 
influencia, que predice la fórmula (23) y de la reobase que predice la (25), 
han podido ser confirmados en experimentos in-vivo, in-vitro e in-silico 
para el caso de estimulación eléctrica con electrodos puntuales en condi- 
ciones en las que la longitud del intervalo de influencia es proporcional a 
la distancia entre el electrodo y la fibra [12] [13] [18]. 

Para estimulación magnética con equipos de tercera generación, al parecer 
todavía no hay datos adecuados. 


(d)La deducción de fórmulas para la curva intensidad-duración correspon- 
dientes a la excitación catódica se basó en la dinámica del umbral para el 
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modo fundamental, considerado como modo dominante y desacoplado de 
los demás modos. 

Cabe preguntarse por el error que se comete al estimar el umbral del modo 
dominante desacoplado en la estimulación catódica. La Fig. 8 muestra los 
resultados en forma de gráfica, para el caso de una estimulación eléctrica 
catódica empleando un electrodo puntual y pulsos rectangulares de co- 
rriente de amplitud /,. 


XXX Simulación numérica con modos acoplados 


——- Aproximación adiabática 


== - Modo dominante desacoplado 


0 10 20 
z (mm) 
Fig. 8. Comparación del umbral para el modo dominante desacoplado con los umbrales 


para dicho modo calculados a partir de una aproximación adiabática y a partir de la diná- 
mica de modos acoplados. 


En abscisas aparece la distancia entre el electrodo y la fibra y en ordenadas 
se encuentra la amplitud umbral. 

Las cruces corresponden a los valores del umbral obtenidos a partir de la 
dinámica con modos acoplados. 

La curva continua, que coincide prácticamente con el resultado de la di- 
námica de modos acoplados, se calculó suponiendo que los modos de or- 
den superior al primero evolucionan en cuasi-equilibrio con este modo 
dominante (aproximación adiabática o por modos esclavizados). 

La curva que corresponde al modo dominante desacoplado coincide con 
las otras dos cuando el electrodo se encuentra a distancias cortas respecto 
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de la fibra, pero se dispone por encima a distancias mayores (sobre-estima 
el umbral). 

Para que los modos de orden superior (3, 5,...) puedan describirse median- 
te una aproximación adiabática la duración del pulso estimulador debe ser 
por lo menos el triple de la mayor constante de tiempo de las que caracte- 
rizan a los modos esclavizados. En este caso dicha constante es la del mo- 
do 3. 

Las amplitudes correspondientes a los modos esclavizados se pueden esti- 
mar mediante las ecuaciones: 


E dB, -AJlt)+ AF, 1, 


2 
1 2 JE 
A 


Sustituidas en la ecuación para el modo dominante acoplado con los demás 
modos, mientras dura el pulso estimulador, se obtiene una ecuación de 


d 
evolución del tipo qa (1)= A a AA e Flo 


A partir de esta ecuación se puede analizar el bloqueo anódico por el cá- 
todo. En ciertas condiciones, cuando la amplitud de la corriente catódica 
continúa aumentando una vez alcanzado el umbral, se llega a una amplitud 
de bloqueo por encima de la cual la fibra ya no resulta excitable. El análi- 
sis puede llevarse a cabo aplicando el mismo método empleado para el 
caso de las fibras no revestidas de mielina [15], [16]. 


4,(1)= 


(e)En esta primera parte no se tuvo en cuenta el efecto de la variable de 
recuperación, que se hace más importante a medida que aumenta la dura- 
ción del pulso estimulador y que permite explicar las otras dos modalida- 
des de excitación, no consideradas en esta parte del trabajo: la estimula- 
ción anódica de apertura y la estimulación catódica de apertura. La aplica- 
ción del análisis modal no lineal al estudio de estos fenómenos es el pro- 
pósito de la segunda parte. 

Finalmente, un comentario sobre el enfoque, fundamentalmente analítico, 
de este trabajo. 

Modelos que no atienden al detalle anatómico y fisiológico pueden ser no 
triviales desde el punto de vista matemático. Se puede trabajar con mode- 
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los bastante realistas utilizando alguno de los potentes paquetes de softwa- 
re disponibles. No obstante, una exploración completa de la totalidad de 
los comportamientos que podrían resultar de interés biomédico, exige se- 
leccionar bien los valores de numerosos parámetros implicados en el mo- 
delo. Esto puede ser muy difícil. En algunos casos, sin información adicio- 
nal, casi siempre se van a omitir combinaciones de parámetros interesan- 
tes. Las fórmulas analíticas como las deducidas en la primera y en la se- 
gunda parte de este trabajo permiten una comprensión más profunda de la 
interrelación entre las variables y los parámetros que aparecen en el mode- 
lo reducido. Este nivel de comprensión no puede obtenerse mediante co- 
rridas de simulación digital realizadas a pura fuerza bruta. Permite diseñar 
mejores y más completos estudios numéricos sobre modelos más comple- 
jos y realistas desde el punto de vista biomédico. 
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neizada. Parte II: Dos variables de estado. 


Roberto Suárez Ántola 
Profesor de Ingeniería y Matemática Aplicada 
Universidad Católica del Uruguay 
8 de Octubre 2801- CP11600, Montevideo, Uruguay 
rsuarez O ucu.edu.uy 


Resumen: Al igual que en la primera parte, en esta segunda parte se emplea una ecuación 
de cable no lineal homogeneizada para describir la interacción entre un campo eléctrico 
externo y una fibra nerviosa mielínica. Aplicando nuevamente el concepto de intervalo de 
influencia del campo sobre el axón, y un modelo simplificado de membrana excitable con 
dos variables de estado, se aproxima la ecuación a derivadas parciales no lineal del mode- 
lo original por un sistema de ecuaciones diferenciales ordinarias en amplitudes de modo. 
Se discuten algunos resultados de simulación digital en el espacio de amplitudes de modo 
y se identifican una amplitud dominante que permite simplificar el orden del sistema de 
ecuaciones diferenciales ordinarias. Se estudia la dinámica del umbral en el sistema de 
orden reducido, empleando los conceptos de conjunto de decaimiento, conjunto amplifi- 
cador y barrera umbral. Se fundamentan las teorías clásicas de dos factores para la excita- 
ción. Se obtienen fórmulas analíticas para las curvas intensidad-duración correspondien- 
tes a la estimulación catódica de cierre y la estimulación anódica de apertura en el caso 
eléctrico, y para las situaciones análogas en el caso magnético. Se discuten las limitacio- 
nes del modelo y se sugieren mejoras. 

Palabras clave: Ecuación del cable homogeneizada, curvas intensidad-duración, umbra- 
les, cronaxia, reobase, parámetro sigma, estimulación eléctrica funcional, estimulación 
magnética, fibra nerviosa mielínica, análisis modal no lineal, estimulación catódica de 
cierre y apertura, estimulación anódica de apertura y cierre. 


Abstract: A nonlinear and homogenized cable equation is employed here, as it was in the 
first part of this work, to describe the interaction between an external electric field and a 
myelin nervous fiber. Applying again the concept of interval of influence of the field over 
the axon, and a simplified model of excitable membrane with two state variables, the 
nonlinear partial differential equation of the original model is approximated by a system 
of nonlinear differential equations in mode amplitudes. Some digital simulation results in 
the space of mode amplitudes are discussed, and dominant mode amplitude is identified 
and used to simplify the system of ordinary differential equations. Threshold dynamics is 
studied in the reduced order system applying the concepts of decaying set, amplifying ser 
and threshold barrier. A foundation is given to the classical two factor theories of excita- 
tion. A mathematical model of strength-duration curves for electric and magnetic stimula- 
tion is constructed. Analytical formulae for strength-duration curves are obtained for 
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cathodic make and anodic brake excitations, in the electric case, and for the analogous 
cases in magnetic stimulation. Some limitations and improvements of the model are dis- 
cussed. 

Key words: Homogenized cable equation, strength-duration curves, thresholds, chronax- 
le, rheobase, sigma parameter, functional electric stimulation, magnetic stimulation, mye- 
lin nervous fiber, non-linear modal analysis, cathodic make and brake excitation, anodic 
make and brake excitation. 


(1)Introducción 


En la primera parte de este trabajo [21] revisamos algunos aspectos de la estimu- 
lación eléctrica y magnética de células excitables con fines terapéuticos e insisti- 
mos en que el mecanismo de excitación o de bloqueo involucra una interacción 
local entre un campo eléctrico externo (ya sea producido por un sistema de elec- 
trodos externos, ya sea inducido en el conductor de volumen de los tejidos por un 
campo magnético variable en el tiempo) y la célula excitable cuya actividad se 
desea estimular o bloquear. 

En el caso de la estimulación eléctrica el campo se produce mediante un arreglo 
de electrodos, ubicados más o menos próximos al tejido cuya actividad se preten- 
de estimular o bloquear. A través de los electrodos se inyecta una corriente eléc- 
trica en el conductor de volumen formado por los tejidos biológicos, con una 
duración y una amplitud que se seleccionan buscando que la estimulación resulte 
exitosa y específica. 

En el caso de la estimulación magnética, la variación temporal del campo magné- 
tico genera un campo eléctrico inducido, y es este último el responsable de la 
estimulación o el bloqueo en la célula blanco. El flujo magnético variable se logra 
inyectando una corriente eléctrica variable en una bobina o conjunto de bobinas. 
La variación del flujo magnético en el conductor de volumen genera una distribu- 
ción espacial de corrientes y campos eléctricos inducidos. 

Como vimos en la primera parte del presente trabajo, el resultado del proceso de 
estimulación eléctrica o magnética en el sistema formado por el dispositivo esti- 
mulador, el conductor de volumen de los tejidos y el elemento o elementos que se 
desea estimular, puede describirse a nivel de caja negra como una corresponden- 
cla entre cada historia de corriente eléctrica aplicada (inyectada en el conductor 
de volumen de los tejidos mediante electrodos o inyectada en las bobinas para 
producir el flujo magnético variable) y una variable binaria A que toma el valor O 
si la estimulación fracasa y 1 si tiene éxito. 

La variable binaria A refleja el carácter de fenómeno de umbral que presenta la 
estimulación eléctrica y magnética. Es posible determinar la correspondencia 


Rev. SCP”— — Volumen18 N*1 pp. 17-46 ISSN 0379-9123 


entre historias de corriente y el valor que adopta Á a partir de medidas eléctricas 
de la respuesta del o los elementos que se desea estimular, u observando manifes- 
taciones externas como contracciones musculares o inhibición de funciones. 

La determinación de la curva intensidad-duración es una herramienta que se utili- 
za en la práctica para llevar a cabo la identificación del umbral que en definimos 
mediante la variable A [15], [20]. 

La medición de los umbrales de estimulación eléctrica o magnética debe hacerse 
en cada caso, porque la curva intensidad-duración depende de la posición y tama- 
ño de los electrodos o las bobinas, de las propiedades del conductor de volumen 
formado por los tejidos biológicos y de la naturaleza y posición de la o las células 
blanco de la estimulación [1], [8], [9], [10], [12], [15], [16], [17], [18], 20]. 

En la primera parte de este trabajo construimos un modelo matemático aplicable a 
la descripción de los fenómenos de umbral en la estimulación eléctrica y magné- 
tica de fibras. El modelo se basa en la hipótesis del intervalo de influencia del 
campo externo sobre la fibra, en una ecuación de cable no lineal homogeneizada 
junto con ecuaciones simplificadas para la unidad de membrana en los nodos de 
Ranvier, y en el método de análisis modal no lineal. 

El concepto de intervalo de influencia permite suponer, solamente hasta que se 
alcanza el umbral, que las variables de estado de la unidad de membrana perma- 
necen en sus valores de reposo en los puntos extremos de dicho intervalo. 

La representación del axón mediante una ecuación del cable homogeneizada fue 
propuesta y fundamentada por Basser con el propósito de describir el comporta- 
miento de la fibra mielínica por debajo del umbral [2]. Como dijimos en la pri- 
mera parte, en este trabajo extendemos su ámbito de aplicación para describir 
perturbaciones de las variables de estado de la membrana hasta alcanzar un esta- 
do umbral. La membrana a la que se hace referencia es una membrana equivalen- 
te, cuyas variables de estado aparecen como variables dependientes del tiempo y 
de la posición en la ecuación homogeneizada. Estas variables de estado emergen- 
tes pueden interpretarse como una especie de promedio de las correspondientes 
variables en los nodos de Ranvier y en los manguitos de mielina. 

El método de análisis modal no lineal permite reemplazar el sistema de ecuacio- 
nes original (que incluye una ecuación a derivadas parciales no lineal), junto con 
las condiciones en los extremos del intervalo de influencia, por un sistema de 
ecuaciones diferenciales ordinarias no lineales en términos de las amplitudes de 
modo correspondientes a las variables de estado de la membrana. Estas amplitu- 
des de modo son funciones del tiempo a determinar, mientras que la dependencia 
del espacio (en el caso de una fibra la variable espacial es la longitud de arco 
tomada a lo largo del axón) viene dada por un conjunto de funciones propias de 
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un operador lineal, y son siempre las mismas. 


El propósito principal de esta segunda parte del trabajo es continuar el análisis del 
modelo matemático construido en la primera parte, teniendo en cuenta ahora el 
efecto de la variable de recuperación sobre la dinámica del umbral. 

A partir de los resultados de algunas corridas de simulación digital y de las consi- 
deraciones teóricas presentadas en [21], confirmamos que el primer modo, que 
había resultado ser dominante para una sola variable de estado (el voltaje a través 
de la membrana), también resulta dominante para el caso en el que se agrega una 
variable de recuperación!. 

Introducimos el conjunto de decaimiento y el conjunto amplificador en el plano 
de amplitudes del modo dominante para el voltaje a través de la membrana y para 
la variable de recuperación, y describimos una barrera umbral que separa, en 
forma más o menos difusa, un conjunto del otro. 

Suministramos una fundamentación a los modelos clásicos de dos factores de la 
excitación eléctrica mediante un electrodo externo y hallamos que en ciertas con- 
diciones cabe esperar que dichos modelos puedan extenderse a la descripción de 
la excitación magnética de fibras nerviosas. 

Hallamos nuevas fórmulas analíticas para las curvas intensidad duración, aplica- 
bles tanto a la estimulación catódica de cierre y anódica de apertura como a las 
situaciones análogas en el caso de la estimulación magnética. 

Comparamos algunas de las predicciones que se desprenden de las fórmulas ma- 
temáticas para los umbrales de estimulación, con los resultados de la determina- 
ción experimental de esos umbrales para diferentes tipos de fibras nerviosas y con 
los resultados de las corridas de simulación digital efectuadas a partir de modelos 
recientes, más complejos y realistas, pero que no permiten el tipo de abordaje 
analítico empleado en este trabajo. 

Finalmente discutimos algunas limitaciones del modelo en su forma actual y su 
extensión para poder tener en cuenta el efecto de las variables de activación sobre 
la dinámica del umbral de estimulación eléctrica y magnética. 


(2) Revisión del modelo matemático de la interacción campo-axón ba- 
sado en una ecuación de cable no lineal homogeneizado 


Y Las corridas de simulación digital, cuyos resultados se presentan y discuten en este 
trabajo, fueron programadas en Matlab y llevadas a cabo por los ingenieros Héctor 
Korenko y Diego Suárez-Bagnasco, empleando en cada caso la selección de parámetros 
sugerida por el autor. 
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Como vimos en [21], mediante un proceso de homogeneización, basado en el 
método asintótico de las escalas múltiples (dos escalas espaciales en este caso, 
una del orden de la longitud del nodo y otra del orden de la distancia entre nodos) 
se puede hallar una ecuación de cable emergente, como una primera aproxima- 
ción, a partir de operar sobre un modelo de la micro-estructura detallada del axón. 
Para los fines del presente trabajo partimos del siguiente sistema: 


Ov ye) 2» 2 Oo par 
(EA la 
a E Os pe 
ya) 
=Llv—y-. 1b 
E (1b) 


En las ecuaciones (1a) y (1b) el campo escalar v(t, s) representa el potencial a 


través de la membrana excitable equivalente (es decir, emergente del proceso de 
homogeneizado), relativo a su valor de reposo, en el instante £ y en la posición 
de abscisa curvilínea sa lo largo del axón (se introduce una abscisa curvilínea 

para poder trabajar con fibras curvadas). El campo escalar w(t,s) representa una 


variable de recuperación de la membrana equivalente. En el cociente E , L re- 


presenta la longitud entre nodos de Ranvier y : la longitud de un nodo. 
Supusimos, y continuaremos suponiendo que el término N(v,w) que representa 


las no linealidades que hacen posible la excitación e incluye el efecto de la varia- 
ble de recuperación w es del tipo: 

N(v,w)=bv*=cv*—a-w (2) 
Los parámetros que aparecen en estas ecuaciones son positivos. 
La constante de espacio 4 del cable homogeneizado se expresa en función de las 


constantes de espacio 4,,,, del manguito de mielina y 4,,del nodo de Ranvier a 


través de la fórmula: A ¡Pe _ E > (3) 

A L Amie L Ano 
La constante de tiempo 7 del cable homogeneizado se expresa en función de las 
constantes de tiempo 7,,, del manguito de mielina y 7,,del nodo de Ranvier 


2 2 
mediante la fórmula: r=[1-£]. d a da “Tno (4) 
L Amic L Ano 


El término que conecta el estado de la membrana equivalente con el campo eléc- 
trico externo, involucra la derivada de la componente del campo externo tangen- 
cial al axón respecto de la abscisa curvilínea del axón. En las condiciones usuales 
de estimulación eléctrica y magnética se puede expresar como un producto de una 
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función del tiempo por una función del espacio (de la abscisa curvilínea s ): 
OE, par(t>5) 
a IN) (5) 
OS 
En la estimulación eléctrica gls) es la corriente inyectada por el electrodo activo 


ig(t), mientras que en la estimulación magnética es la derivada respecto del tiem- 


po de la corriente en la bobina e (1). 
dt 


En ambos casos F(s) representa la variación espacial de la componente del cam- 


po externo paralela al axón, por unidad de longitud tomada a lo largo de la fibra. 
Lo denominaremos factor de forma. 

Por lo general, el proceso de estimulación, desde el reposo y hasta que se alcan- 
za el umbral es un fenómeno espacialmente localizado: fuera de una región de 
influencia acotada es posible despreciar la perturbación que el campo externo 
produce en las variables de estado de la membrana y asumir que estas últimas 
permanecen fijas en sus valores de reposo. Esto sugiere introducir el intervalo de 
influencia, admitiendo que F(s) se anula fuera de un intervalo acotado de longi- 


tud /. Si admitimos esta hipótesis como una primera aproximación, podemos 
suponer que los campos v(+,s) y w(r,s) se anulan en los extremos del intervalo de 


influencia. Si el origen se pone en el punto medio de dicho intervalo: 143) =0 


y u(13) =0 Como vimos en [21], el próximo paso se basa en las propiedades 


del operador L|..]=-? 


vlt, 5). 

Se obtienen los siguientes valores propios y funciones propias de ese operador 

(que forman un conjunto completo de funciones normalizadas y ortogonales): 
2,2 


o == siendo / la longitud del intervalo de influencia y n=1, 2, 3,...; 


que aparece en la ecuación (la) actuando sobre 


Oi 
Os? 


o, (s)= 2 cos) cuando n es impar y q,(s)= (sat) cuando es par. 
" / E y / 


A continuación planteamos desarrollos aproximados en términos de las funciones 
propias dependientes de la posición sobre el cable, con amplitudes de modo de- 
pendientes del tiempo: 
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5lt.s)> Y a,(0)-0,(s) (6 a) 0) Za ldol) 6D) 


Sustituyendo (6 a) y (6 b) en (la) y (1b) y proyectando sobre las primeras M fun- 
ciones propias e igualando a cero esas proyecciones, se obtiene un sistema de 
ecuaciones diferenciales ordinarias no lineales para las amplitudes de modo: 


dA 212 4? AU AS 4 2 
2 = (+ nó) A, +00 0... AJA. -Le Y 0), A, ALA, La B, +40 E, gt) 
dt L L asi pq p.q L dare pq Pp q L 
(7 a) 
¿ee =2.(a, -y B,) (n=1, Ze Bor M) ( b) 
o L ñ 
2 
En estas expresiones se define Oapy= Son ts)-Pp(s)-p,(s)- ds , mientras que 
14 


2 


2 
OGapar= ) pu [s)-p,(s)- p,(s)- p,(s) ds Estos coeficientes dependen solamente de las fun- 


ciones propias del operador: son independientes de la naturaleza del sistema for- 
mado por el dispositivo estimulador, el conductor de volumen y la fibra blanco de 
la excitación. Sus expresiones analíticas pueden obtenerse mediante el uso de 
paquetes de software para cálculo simbólico o recurriendo a tablas de integrales. 

Las proyecciones del factor para la estimulación eléctrica Oo magnética vienen 


2 
dadas por: F, = [on(s) Fl): ds (n=1, 2, 3,..., M) (8) 
+ 


En la referencias [11], [15] y [16] se pueden hallar los valores numéricos de los 


coeficientes O,» Onpgr Y también, para un caso particular, los valores de los 


proyecciones F,. En todos los casos estas últimas proyecciones tienden a cero 


cuando el índice n del modo aumenta. 

El sistema de ecuaciones original, que incluye una ecuación a derivadas parciales 
no lineales y describe la dinámica de la membrana del cable homogeneizado, 
queda así sustituido por un sistema de ecuaciones diferenciales ordinarias para las 
primeras M amplitudes de modo A, (+) y B,(+) correspondientes a las expresiones 


aproximadas (fórmulas (6 a) y (6 b)) de las variables de estado. Cada ecuación de 
evolución correspondiente a una amplitud de modo presenta un término forzante 
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proporcional a la corriente ¡¿(+) en el electrodo en el caso de la estimulación eléc- 


trica, o a la velocidad de variación de la corriente en la bobina te (+) en el caso 
dt 


de la estimulación magnética. 
En el factor de proporcionalidad aparece la proyección F, del factor de forma 


dada por la ecuación (8). Si esta proyección es nula o al menos despreciable, el 
modo solo puede activarse a través de su acoplamiento con los otros modos, al- 
gunos de los cuales se verán forzados por unas proyecciones no nulas. 


(3) Análisis modal no lineal de la ecuación del cable homogeneizada: 
dinámica del umbral para una membrana con dos variables de estado. 


(3.1) Simulación digital de la estimulación utilizando tres amplitudes de mo- 
do para cada variable de estado de la membrana excitable. 

Comencemos revisando la Fig. 1, ya presentada en la primera parte de este traba- 
jo. Muestra los resultados de una simulación digital efectuada a partir de las 
ecuaciones (10 a) para M=5, poniendo «a =0 de modo de eliminar el efecto de la 
variable de recuperación. 


1 


Amplitud de modo 


% 0.5 1 


tiempo (ms) tiempo (ms) 
Fig. 1: Dinámica forzada por un pulso catódico (parte izquierda) y dinámica libre una vez 
finalizado el pulso estimulador (parte derecha) con tres amplitudes de modo [|]. 


La simulación se llevó a cabo para el factor de forma simétrico de un electrodo 
puntual, que activa solo tres amplitudes de modo (A;¡, 43 y 45). La dinámica 


de la fibra se describe entonces en un espacio de estados de tres dimensiones. 


Para hacer la simulación digital se tomó: “b=1 4.-1 r=3ms A4=1mm y una 
L L 8 


distancia electrodo fibra de z=20mm(dato necesario para poder calcular las 
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proyecciones F, F3 Fs del correspondiente factor de forma []. Además se 
necesitan los valores de los parámetros 6G,¡¡=1.2, 6,¡¡¡=1.5que dependen sola- 


mente de la primera función propia del operador espacial, y, (s)= pa se 22) 
Como g(+) se eligió un pulso rectangular de amplitud negativa, que en este caso 


en particular, en el que la estimulación es eléctrica, corresponde a un pulso cató- 
dico. La duración del pulso se fijó en 1 ms, de modo que resultara corto respecto 
de la constante de tiempo de la membrana de la fibra homogeneizada. 

La figura de la izquierda muestra la activación de las amplitudes de los tres mo- 
dos de índice impar a partir de sus valores de reposo, mientras dura el pulso esti- 
mulador. 

La figura de la derecha muestra la evolución de esas mismas amplitudes de modo 
una vez que el pulso estimulador ha finalizado. 


La Figura 2 muestra el resultado de la simulación digital con el mismo pulso apli- 
cado a través del mismo cátodo puntual a la misma distancia de la fibra, con los 


mismos parámetros de la membrana homogeneizada, excepto que ahora qua = 1, 


T . o ; : 
2. E=01 y y=0.1 porque la variable de recuperación se ha incluido en el es- 
w 
tudio. 
5 4 
€ 
y 3 
1 
5 3 3 
E o 
3 3 
3 [7] 3 p 
És esa 
32 p 25 
3 > Lo] S 
En E > A 
o E 
dy A 


0 6 12 0 6 12 
Tiempo adimensional o = - Tiempo adimensional o = z 

Figura 2: Dinámica libre con seis amplitudes de modo (tres para el campo de voltaje a 
través de la membrana, a la izquierda, y tres para el campo de la variable de recuperación, 


a la derecha) luego de finalizado el mismo pulso catódico utilizado para la Figura 2. 
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Los valores asintóticos de las tres amplitudes de modo que aparecen en la parte 
derecha de la Figura 1 corresponden a un estado excitado de la fibra, espacial- 
mente no uniforme y estable, debido a que la variable de recuperación se mantie- 
ne congelada en su valor de reposo. La amplitud A, del modo fundamental es 


dominante respecto de las otras dos amplitudes de modo, pese a que al final del 
pulso estimulador su amplitud es inferior a las amplitudes de los otros dos modos. 
La amplitud del pulso estimulador es lo bastante grande como para que, una vez 
finalizado, A, tienda asintóticamente a un valor no nulo, arrastrando consigo a 


las otras amplitudes, que se mantienen en cuasi-equilibrio con la amplitud de 
modo dominante. 

La representación de las seis amplitudes de modo que aparecen en la Figura 2 
comienza justo cuando el pulso estimulador ha finalizado. Después de una excur- 
sión más o menos extensa (cuya etapa final no se ha representado en la figura), 
todas las amplitudes de modo terminan por regresar al estado de reposo, debido al 
efecto de la variable de recuperación. Pero tanto A, como B,son dominantes 
respecto de las demás amplitudes de modo correspondientes a la misma variable 
de estado de la unidad de membrana. 

Estos y otros resultados análogos de simulación digital? sugieren que el primer 
modo resulta ser dominante, también cuando la variable de recuperación de la 
membrana se incluye en la dinámica. Admitiendo que esto es así, en una primera 
aproximación se lo puede desacoplar de las otras amplitudes de modo. 


(3.2) Simulación digital de la dinámica en el plano de estados formado 
por las amplitudes de modo dominantes. 


Con esta simplificación, a partir de las ecuaciones (7 a) y (7 b) obtenemos una 
dinámica no lineal en un espacio de estados bidimensional que en principio po- 
dría aplicarse tanto a la estimulación eléctrica como magnética, si el intervalo de 
influencia y la proyección del factor de forma se interpretan adecuadamente: 


dA PR 

n= (er? A +E-b:0, APEC: Aja B +2 E 80) (9 a) 
dB, _u 

Tr, =L.(A -y-B, (9 b) 
w di L ( 1 Y ,) 


El tipo de sistema de ecuaciones diferenciales ordinarias, al cual pertenece nues- 
tro sistema (9 a)-(9 b), se pueden denominar de Bonhófer-van der Pol, en honor a 


? Las corridas de simulación digital de (7 a) y (7 b) para el caso de una fibra desprovista 
de revestimiento de mielina, excitada a través de un electrodo puntual y teniendo en cuen- 
ta 17 modos, muestran la dominancia del modo fundamental sobre los demás modos [11]. 
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los investigadores que en la segunda y tercera décadas del siglo pasado constru- 
yeron casos particulares y los aplicaron al estudio de problemas de excitabilidad. 
En [7] puede hallarse una investigación de las propiedades de otro caso de este 
tipo de sistema, empleado como modelo matemático simplificado para describir 
la dinámica de una porción de membrana excitable uniformemente polarizada. 
En el caso estudiado por Fitz-Hugh en [7], las variables de estado son el voltaje a 
través de la membrana (en lugar de nuestra A, ) y una variable de recuperación (en 
lugar de nuestra B, ), mientras que los parámetros son 7,,,7,,, b, C, AL y y, O sea los 
parámetros de lo que se suele denominar una unidad de membrana. 
Los parámetros adicionales GQ, ¡¡,Q, 111, 4,4, 11 y Lque aparecen en nuestras ecua- 
ciones () y (), en el modelo estudiado por Fitz-Hugh no aparecen. 
Los resultados de simulación digital que se presentan a continuación representan 
el comportamiento de la fibra abandonada a sí misma (dinámica libre), una vez 
alcanzado un cierto estado inicial: supondremos entonces que g(t) se anula a 
partir del instante inicial. El estado inicial del sistema puede ser cualquiera de los 
accesibles por estimulación mediante electrodos o bobinas, o simplemente puede 
haber sido establecido como tal en forma arbitraria, con el fin de estudiar la di- 
námica en el plano de fases (A, B,). Para la simulación digital es conveniente 
eliminar una de las escalas de tiempo introduciendo un tiempo a-dimensional: 
í : SS A 
o == Como la constante de tiempo 1 cabe presumir que se mantiene incambiada 
7 
durante la estimulación, la utilizamos como referencia de tiempo básica. En este 
caso, para una dinámica libre las (9 a) y (9 b) se reducen a: 


dA RP 
> = -(1+ 7? DAMA Ec Al Lab (10 a) 
dB; T 2) 

= . -(A, -B 10b 
> (E) (4, =7-B,) (10b) 


El parámetro a-dimensionado £ (5) es siempre menor que 1 debido a que la 


Tay 


escala en la que se activan los canales que producen la recuperación de la mem- 
brana es mayor que la escala correspondiente a la constante de tiempo de la 
membrana. 

Cuanto menor sea este número adimensional, mayor va a ser la rigidez del siste- 
ma de ecuaciones diferenciales no lineales que describen la dinámica del umbral 
en términos de las amplitudes dominantes. 

Teniendo todo esto en cuenta quedan aún por fijar valores para siete parámetros: 
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ES 
( 


Y 17 17 
—e. 2.b —C 
pe L L 


He 
S 
— 


Fijamos los siguientes valores de = 0.00833 a «b=1 E :-c=0.125 = -q=1 
/ 


y=0.1 
En la de la parte izquierda de la Figura 3 se pueden ver algunas órbitas represen- 
tativas para esta primera elección de los parámetros. En las corridas de simula- 


ción correspondientes se eligió dE Pp==:—= 


En la parte derecha de la Figura 3 todos los parámetros, con excepción de f son 
iguales a los de la parte izquierda. Pero ahora hay una diferencia significativa 
entre la escala de tiempo de las variables de recuperación y la escala de tiempo de 


17, 


la membrana debido a que se tomó L - f' igual a 0.1. Las órbitas representativas 


muestran que se produjo una modificación fundamental en la dinámica del um- 
bral. 


PLANO DE FASES (A1,B1), alfa=1, beta=0.3, gamma=0.1 
r Y Y —r 7 Y 7 Y 


a o 
SS 


- A E 5 6 s 2 4 o 1 2 3 4 BL 8 
A1 
Figura 3. Órbitas en el plano de fases para el primer juego de parámetros (izquierda) y 
para el segundo juego (derecha). 


PLANO DE FASES (A1,B1), alfa=1, beta=0.1, gamma=0.1 


En ambas figuras se han representado las dos isóclinas del sistema dinámi- 


co. La isóclina correspondiente a B, _ y, es la recta B,= 1 A,» Independien- 
dt Y 


te del cociente 4 entre la longitud del nodo de Ranvier y la distancia inter- 
E 
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nodal. La isóclina correspondiente a a = (, es una curva sigmoidea hori- 
t 
. L 2 R 2 3 
zontal dada por: B,=-——-(l+x JA, +b-8,,,A¡ —c:8 111: Aj (11) 
au £ 

En líneas generales, las órbitas que pasan por o parten de puntos próximos a la 
rama media de la curva sigmoide, a uno y otro lado de la misma, presentan dos 
comportamientos diferentes. 

Los puntos a la izquierda de la rama media originan órbitas que se dirigen a la 
rama izquierda de la sigmoide y descienden muy próximas a ella hacia el origen. 
Cuando se examinan las funciones (A (2), B (+) que corresponden a esas órbitas, 


se advierte que no presentan un comportamiento asignable a potenciales de ac- 
ción (PA). Las variables retornan al estado de reposo (el origen del plano de las 
fases determinado por las variables (A, B,)) sin presentar la excursión caracterís- 


tica de los PA. Si estas trayectorias representaran el resultado de un intento de 
excitación de la fibra, los pulsos elegidos habrían sido sub-umbrales y la excita- 
ción no se habría producido. 

Por su parte, los puntos que se encuentran a la derecha de la rama media de la 
sigmoide originan órbitas que se dirigen primero hacia la rama derecha de la sig- 
moide, ascienden cerca de ella y luego la abandonan dirigiéndose hacia la rama 
izquierda de la sigmoide para retornar al origen descendiendo próximas a ella. 
Cuando se examinan las funciones (a (0), B (2) que corresponden a esas órbitas, 


se constata que presentan las características asignables a los PA en el marco del 
presente modelo. Ahora, cuando se dirigen y alcanzan la rama derecha de la curva 
sigmoide, para luego regresar hacia el estado de reposo como se describió pre- 
viamente, las variables presentan la excursión característica de los PA. Si estas 
trayectorias representaran un intentos de excitación eléctrica o magnética de la 
fibra, los pulsos elegidos habrían resultaron supra-umbrales y la excitación (géne- 
sis de una respuesta en forma de PA potencial de acción) se habría producido. 

Las corridas de simulación permiten identificar un conjunto de órbitas de transi- 
ción entre las que presentan las características de un PA y las que claramente no 
las presentan. Este conjunto de transición es particularmente evidente en la parte 
izquierda de la Figura 3. 


(3.2) Dinámica del umbral para el modo dominante desacoplado: con- 
junto de decaimiento, conjunto amplificador, curva umbral y regiones 
accesibles mediante un sistema de estimulación. 


Estos dos tipos de comportamiento de las órbitas en el plano de fases, cuya dife- 
rencia es mucho más nítida en el caso que muestra la parte derecha de la figura 3, 
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que se aproxima al fisiológico y en consecuencia presenta un mayor interés para 
este trabajo, sugiere introducir una curva umbral (pseudo-separatriz umbral) 
como se ve en la parte izquierda de la figura 4. Dicha curva umbral que delimita 


dos regiones: el conjunto de decaimiento y el conjunto amplificador. 
A E 


B1 


2 dA1/dt = 0 
Conjunto de 
decaimiento 


Conjunto 
amplificador 


Pseudo-separatriz umbri 


A1 - 2 4 0 1 2 3 4 E 6 
Figura 4. A la izquierda se representan las isóclinas, los conjuntos de decaimiento y am- 
plificador, la pseudo-separatriz umbral y las regiones del plano de fases accesibles me- 
diante la estimulación eléctrica o magnética. A la izquierda se muestra la barrera umbral. 


Si un pulso g(t) (de corriente ¡,(t) en la estimulación eléctrica y de velocidad de 


variación de la corriente en la bobina E (+) en el caso de la estimulación magné- 
dt 


tica) lleva el estado de la fibra desde el reposo (el origen del plano de las fases) 
hasta un estado situado en el conjunto de decaimiento y en ese momento cesa el 
pulso, el estado evoluciona en el conjunto de decaimiento retornando al origen, y 
no se produce un PA en la fibra: la excitación fracasa. Pero si un pulso de corrien- 
te lleva al estado de la fibra desde el origen hasta cruzar la curva umbral para 
entrar al conjunto amplificador, y luego cesa el pulso, se desencadena un compor- 
tamiento del estado del sistema que corresponde a la aparición de un PA: la esti- 
mulación resulta exitosa. 

En la parte izquierda de la Figura 4, además de la curva umbral y los conjuntos de 
decaimiento y amplificador, se han indicado con regiones sombreadas los conjun- 
tos accesibles desde el origen (estado de reposo de la fibra), mediante pulsos de 
corriente aplicados a través de un sistema de estimulación eléctrica o magnética. 
El sector angular sombreado que aparece a la derecha y por encima del eje Aj, es 
el accesible mediante pulsos catódicos en el caso eléctrico o sus equivalentes en 
el caso magnético. El otro sector sombreado es el accesible mediante pulsos 
anódicos o sus equivalentes en la estimulación magnética. 

Para ver esto, supongamos que la membrana de la fibra parte del reposo (es decir, 
del origen del plano Aj, B;). De la (9 a) se desprende que para un pulso catódico 
A, abandona el origen en sentido del eje A, positivo. En el caso de un pulso 
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anódico A, abandona el origen en sentido del eje A; negativo. De la (9 b) se des- 
prende que en el caso de un pulso catódico B, abandona el origen siguiendo a Aj, 
manteniéndose siempre A,>y-B,. Entonces el estado del sistema se mantiene 


siempre por encima del semieje A, positivo y por debajo de la recta p, se A» €S 
7 


decir en el área sombreada que se encuentra en el primer cuadrante del plano de 
las fases. Empleando un razonamiento análogo, se ve que para un pulso anódico 
el estado del sistema queda en la zona sombreada del tercer cuadrante de la mis- 
ma figura. Como consecuencia solo se pueden alcanzar, mediante sistemas de 
estimulación, partiendo del estado de reposo de la fibra, los puntos de la curva 
umbral que pertenecen a uno de los dos conjuntos accesibles. Esos puntos de la 
curva umbral son los estados umbral de la fibra que pueden ser alcanzados me- 
diante el sistema de estimulación considerado. 

Finalmente, en la parte derecha de la Figura 4 aparecen la curva sigmoide y la 
recta tangente a dicha curva, en el punto en el cual la rama media de la sigmoide 


corta el eje A,. A partir de la ecuación (11) que relaciona B,con A, se pueden 


b-6, 
de An A 


hallar los puntos de intersección con el eje 4;: 
A =0 Ar,= 4 
209111 


E b-0,11 
e 20491111 


A 2 e 
En estas fórmulas , (2:45 Iii? 4| debe ser positivo para 
2001111 ue La 
L Gu 


que existan los últimos dos puntos de intersección. 
La amplitud A,,, (0) donde la rama media interseca el eje de abscisas correspon- 


de a la amplitud umbral para la estimulación hallada en la primera parte de este 
trabajo [21]. La barrera umbral viene dada por la fórmula: 


a BO (a 4,60) (12) 


1 
1 


Esta recta, denominada barrera umbral, es útil para distinguir en forma apro- 
ximada las órbitas sub-umbrales (que corresponden a una excitación fracasada) 
de las órbitas supra-umbrales (que corresponden a una excitación exitosa): las 
órbitas sub-umbrales se mantienen siempre en el semiplano que se encuentra por 
encima de esa barrera umbral (por tanto no la intersecan), las órbitas supra- 
umbrales siempre intersecan la barrera umbral y pasan al semiplano situado por 


* La denominación barrera umbral se debe a Richard Fitz-Hugh. La introdujo en su análi- 
sis de la génesis de los PA en una unidad de membrana uniformemente polarizada [7]. 
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debajo de la barrera, mientras que la órbita que podemos denominar justo- 
umbrales es tangente a la barrera umbral. 


(4)Deducción de fórmulas analíticas para las curvas intensidad dura- 
ción en condiciones de excitación catódica de cierre y de excitación 
anódica de apertura. 


La posibilidad de introducir las barreras umbral en los planos de las amplitudes 
de modo, para distinguir órbitas sub-umbrales, justo umbrales y supra-umbrales, 
permite generalizar el método empleado en la primera parte de este trabajo, que 
describía el proceso de polarización de la fibra mediante una sola variable de 
estado ( A,) al caso en el que se emplean dos variables (A, y B,) para describir 


esa polarización. 
Como vimos en [21], comenzamos linealizando la ecuación de evolución de la 


amplitud de modo A, en torno al origen. Suponemos que el sistema estimulador 


comienza a polarizar la fibra a partir del estado de reposo de esta última, es decir 
desde A,=0. Introdujimos luego la amplitud umbral A, , y supusimos que la di- 


námica linealizada se podía aplicar para describir la polarización de la fibra hasta 
que se alcanzaba esa amplitud umbral. Una vez alcanzado o apenas superado el 
umbral, se produce el PA y la dinámica linealizada resulta inaplicable. 

Ahora, con dos variables de estado, también comenzamos aproximando lineal- 
mente la dinámica por encima de la barrera umbral: 


dA 2 
A, 0 B +2 E 80) (13 a) 
dB, _4 
=£.(A -—v-B 13 b 
Tr eL (4, Y 5) ( ) 


En lugar de introducir una amplitud umbral, introducimos una curva umbral en el 
plano de fases. Para construir esta curva, consideramos el tramo de la única órbi- 
ta que, viniendo desde la izquierda del plano de las fases llega tangente barrera 
umbral. Ese tramo de órbita se puede hallar en forma aproximada resolviendo el 
sistema linealizado (13 a) y (13 b) poniendo g(t) =0 de modo de trabajar con la 


dinámica libre de las amplitudes de modo. Con este tramo de órbita y la semirrec- 
ta de la barrera umbral que nace del punto de tangencia y se dirige hacia arriba y 
a la derecha, formamos una aproximación a la curva umbral. 

Una vez construida una aproximación a la curva umbral como la descripta, el 
plano de fases queda dividido en dos regiones. 

Los puntos de la región situada por encima forman el conjunto de decaimiento. 
Los puntos de la región que se encuentra debajo de la curva umbral, forman el 
conjunto amplificador. 
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En el croquis que aparece en la Figura 5, la línea gruesa marcada con doble flecha 
representa la curva umbral formada por un tramo de órbita y la mencionada semi- 
rrecta que forma parte de la barrera umbral. El conjunto de decaimiento se en- 
cuentra por encima y el conjunto amplificador se halla debajo de la curva umbral. 


Órbita forzada 
por el sistema 
de estimulación 


Inicio del potencial 
de acción 


| 


. 
Órbita del sistema As =0 
librado a si mismo 
Figura 5. Aproximación a la curva umbral (curva con doble flecha), conjunto de decai- 
miento (área punteada y sombreada) y conjunto amplificador (área en blanco). Se ha 


representado un ejemplo de estimulación anódica de apertura. 


En el ejemplo representado en el croquis de la Figura 5, el sistema estimulador 
(eléctrico o magnético) lleva el estado de la fibra, a través de la región accesible a 
la estimulación anódica, desde el estado de reposo hasta un punto del conjunto 
amplificador, que se asume muy próximo a la curva umbral. La dinámica apro- 
ximada viene regida entonces por las ecuaciones (13 a) y (13 b). A partir del 
momento en que cesa la estimulación, comienza una órbita regida por la dinámica 
libre. El PA se inicia cuando la órbita alcanza la barrera umbral. 

La estimulación catódica produce órbitas que parten del origen, atravesando la 
región accesible a ese tipo de estimulación, y se dirige hacia la parte recta de la 
curva umbral, situada por encima del eje de abscisas. 

Cuanto más intenso es el pulso estimulador, más próxima al eje de abscisas se 
ubica la órbita forzada por el sistema estimulador, eléctrico o magnético, ya sea 
moviéndose hacia la izquierda por la región accesible a la estimulación anódica, 
ya sea moviéndose hacia la derecha por la región accesible a la estimulación ca- 
tódica. 

Para hallar las curvas intensidad-duración que permiten caracterizar el umbral 
para la génesis de PA, debemos en primer lugar integrar las ecuaciones (13 a) y 
(13 b). 
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Para ello introducimos un vector de estado x(1) y las llevamos a la forma matri- 


cial siguiente: o x(1)=A-x(0)+ ele)-b (14) 
t 


2 ] == é 4 
A, (15) ¿Eva sepa (16) Ar A 2 ] a7) 


0 A= 


>=! 
1! 


B, 


Puesto que la fibra parte del reposo, x(0)=0 Aplicamos un pulso g(r) durante el 
intervalo [0,17]. La solución de en estas condiciones es: 


02). l glu)-(e(=0% ¿du (18) 


Si suponemos que describe en todo instante la dinámica del sistema (lo cual 1m- 
plica mantenerse en el semiplano de decaimiento) esta fórmula se verifica para 
todo t positivo, inclusive mayor que la duración 1, del pulso estimulador. 


Para que una órbita desencadene un PA, el vector de estado x(t) debe cruzar 


(proveniente del conjunto de decaimiento y forzado por el sistema de estimula- 
ción) la recta umbral, o al menos ser tangente a la misma. 

Mediante un pulso catódico se llega al umbral cuando la órbita interseca la recta 
umbral. Si justo cuando 1=1, el estado x(+) del sistema coincide con un punto 


de la recta umbral, el pulso de corriente se considera un pulso umbral. 
El pulso es supra-umbral si esta coincidencia se produce antes de t=1p 


El pulso catódico es sub-umbral si el estado del sistema nunca coincide con un 
punto de la recta umbral. 

Mediante un pulso anódico, en el marco de este modelo que emplea solo el modo 
dominante, no se puede llegar al umbral durante el tiempo en el cual ese pulso 
está activo. El comportamiento de la órbita del sistema abandonado a sí mismo 
(después de finalizado el pulso anódico, t>t,) determina si se alcanza o no el 


umbral. Un pulso anódico es justo umbral cuando la órbita que sigue el sistema 
abandonado a sí mismo es tangente a la recta umbral. 

Si la órbita interseca a la recta umbral dirigiéndose hacia el interior del conjunto 
amplificador, el pulso anódico es supra-umbral. Finalmente, es sub-umbral si la 
órbita del sistema abandonado a sí mismo no alcanza la recta umbral. 

Entonces, la condición necesaria y suficiente para que en algún instante el estado 
x(1) pertenezca a la barrera umbral es ¡"x(t)= p , donde p es la distancia entre 


el origen y la barrera umbral, ¡ es un vector unitario columna en el plano de 
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estados que es ortogonal a la barrera umbral y apunta hacia el interior del conjun- 
to amplificador, y 1” es el vector unitario traspuesto (versor fila). Ver Figura 6. 


:d n 
JO Barrera umbral 


4 
-Aá-PpA+=+>—>— >> ¿Barrera umbral 
Figura 6. A la izquierda se puede ver la dinámica libre linealizada en el conjunto de de- 
caimiento, la barrera umbral y la curva umbral aproximada. A la derecha se muestra la 
construcción geométrica que se utiliza para deducir la condición que debe cumplir un 
pulso justo umbral producido por un sistema estimulador, eléctrico o magnético. 


Calculando YA (0) = ] a partir de (11), las componentes según los ejes de 


abscisas n, y de ordenadas n, del versor 1 y la distancia p entre el estado de repo- 


so y la barrera umbral vienen dadas por las fórmulas: 
da, (0) 
103) dB, (19 b) p= 


E) 


Teniendo en cuenta lo anterior un pulso estimulador eléctrico o magnético será un 


AO), 


pulso umbral siempre y cuando: max hiT5(0)= P Q1) 


te|0, +00 
Esto asegura que el estado no abandona el conjunto de decaimiento hasta que 
eventualmente su órbita llega a un punto de tangencia a la barrera umbral. De ahí 
en más la dinámica no se rige por el sistema (14) y de acuerdo con el modelo, se 
genera un PA. 
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2 
Sustituyendo (18) en (19) y dividiendo ambos miembros por ar ne la 
T 


condición para que un pulso sea umbral se puede re-escribir así: 


14 
max fl G(u): glu): as =Cuo (22) 
te[O,+00) 0 
ner Pp 
En esta expresión: G(t)==_" (23) Cuo=7 === — (4) 
ne E A ] -T- 
+ . F a n e 
T 
A partir de las fórmulas (16) y (19 a) se deduce que ¡”¿, = _—— (5) 
| dB, 
En consecuencia, de (20) y (25) se desprende: MES Ars (0) (6) 
U,O 2 


La fórmula (23) puede re-escribirse de esta forma, donde — A, ,— 4, son los valores 


propios de la matriz A : G(t)= rn ¿A Z 7 te (27) 


La figura 7 muestra un ejemplo con valores propios reales [19], [20]. 


LO 


GÍt) 


[| 
tu tm tiempo 


Figura 7. Gráfico de la función de respuesta al impulso donde aparecen el tiempo útil 


hasta la reobase f, y el instante f,, donde la función presenta su mínimo. Esos instantes 


jn2a y t, =2L, 


4-4) A 


se pueden estimar a partir de las fórmulas: *, = ( 


Rev. SCP”— — Volumen18 N*1 pp. 17-46 ISSN 0379-9123 


Calculamos los valores propios a partir de la definición (fórmula (17)) de A, 


2 2 
definiendo los parámetros £= A A NE E? -4(4) O 
T 1580 E Ej 7%, 


w 


1 1 
Los valores propios — A, = > le JA ) y =A,= > le + YA ) pueden ser reales, 
en cuyo caso 0<2A<A,, o bien pueden se  complejo-conjugados 


2=3(Exi-0) A=3(6-1:0) con 0=xV—A cuando A<0 


2 
En este último caso: GÍ(s) =-—,/l + (E ja senío «f— $) siendo tg p= q 
0) a 


1 
A partir de la Figura 7 se ve que para un pulso catódico J Gl: —u): 8. (u)- du 


aumenta mientras el pulso está activo y 1<f,. Una vez que el pulso cesa, co- 


mienza a decrecer. Como consecuencia, para un pulso catódico de duración 
lp 

1, <1, la condición de ser umbral se reduce a fal, -u): 2g.(u)- du=C,.0 
0 

La función GÍls) es a-dimensionada. Depende de los parámetros del modelo de 


fibra excitable homogeneizado y de la longitud del intervalo de influencia (a tra- 
vés de este último, depende del sistema de estimulación). 


Si el pulso es rectangular, su amplitud umbral [e. | se relaciona con su duración 


a través de la fórmula, que podría ser aplicada por igual a la estimulación eléctri- 
ca y a la magnética: 


[e] = Lo (1,<1,) Q8a)  [g.], = =[g.L,. (tp >1,) (28D) 
[0u): cu E lo du 


0 
En (28 b) [e.l, pes el valor reobásico (mínimo) del umbral que se alcanza para 
una duración del pulso que supere al tiempo útil hasta la reobase f, . 

De (28) se desprende que lim,, ol tp-L8. |= Cuo 

En el caso de la estimulación eléctrica [e. | debe interpretarse como la amplitud 


umbral de un pulso rectangular de corriente inyectado a través del cátodo y C,, , 


es la denominada carga umbral límite [3], [5], [7]. En el caso de la estimulación 
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magnética, [e. Ml debe interpretarse como pendiente umbral de una rampa lineal 
de corriente inyectada en la bobina estimuladora mientras que C,, ¿debe interpre- 


tarse como corriente umbral límite. La constante de tiempo para la estimulación 


[e] (c,) 


Cu. (29) y el parámetro sigma O. === (30) también 
8.1. , [e.l.., 
quedan determinados en función de los parámetros del modelo. 
Un análisis de la excitación anódica, análogo al realizado en la referencia [15] 
utilizando la funcional de excitación para reformular las teorías fenomenológicas 


clásicas de la excitación, conduce a la siguiente expresión para la relación entre la 
C 
intensidad umbral ES li y la duración del pulso: le. ], (t.)e => — 61) 


catódica f,.= 


Aquí 1, se determina hallando la menor de las raíces de la ecuación 
GÍ(t,)= GÍr, 4) En el caso que muestra la figura 7, hay una sola raíz, pero 
cuando los valores propios son complejos GÍ(s) oscila con amplitud decreciente y 
hay más de una raíz. También para el caso representado en la Figura 7, la reoba- 
se anódica ES Jl p Se relaciona con la reobase catódica le. IN r A través de la rela- 

j GÍt)- dt 
ción: lec le La > 

fate). ae 


ta 


(32) 


Las correspondientes fórmulas analíticas para la constante de tiempo f_, y el pa- 


rámetro sigma O, para la estimulación anódica se pueden hallar en forma similar 


a la utilizada para la estimulación catódica. 


(5) Fundamentación de las teorías clásicas de dos factores para la es- 
timulación eléctrica y extensión a la estimulación magnética. 


En esta parte nos proponemos re-escribir las ecuaciones linealizadas (13) 
para obtener: (1) Una versión fundamentada de las denominadas teorías de 
dos factores para la excitación eléctrica mediante electrodos externos y (2) 
Una extensión del formalismo de esas teorías, en su formulación original 
meramente fenomenológicas, de modo que la nueva versión que construi- 
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mos resulte aplicable a la estimulación magnética. 
Para comenzar, introducimos las nuevas variables de estado V y U a través de 
las fórmulas siguientes: 


dA, ¿L0) 
V-V,=4, (33 a) U-U,=-=— -B, (33 b) 
dB, 
nde A 1 sá 
Definimos los nuevos parámetros: k == l+x 7 (34 a) 
7 
17) 
72 da, (0) 
o E (34 b) yA 0, (34 c) 
Ki dA,,(0 a dB, 
dB, 
2 7? 
tE) Boy (34d) =p (34 e) 
7 


Teniendo en cuenta lo anterior, las ecuaciones (13a) y (13 b), que descri- 
ben la evolución forzada del estado de la fibra en el conjunto de decai- 
miento, adoptan la siguiente forma: 


dV 
a ki (V-V¿)+k,,(U—=U¿)+ K -2(t) G5 a) 
dU 
de k21 (V-Va)+k,,(U -U,) (35 b) 
Por su parte la barrera umbral, representada por la ecuación (12), se puede re- 
escribir así: e (B,)= An (0) 24:40) B, Entonces el estado (A (0), B (1), 
1 
variando en el conjunto de decaimiento, alcanza un estado umbral cuando 
A() =A,, (0)+ 0 -B, (s) Si igualamos U =A,, (B,) U,=A,., (0) 
1 


V, =0 y tenemos en cuenta las definiciones (33 a) y (33 b), se alcanza el umbral 
de disparo de un PA cuando las nuevas variables verifican la igualdad: 

V=U (36) 
Cuando se interpretan en el marco de la estimulación eléctrica, las ecuaciones 
(35 a), (35 b) y (36) no son otra cosa que una versión muy general de las teorías 
fenomenológicas clásicas de dos factores [3], [4], [7], [8], [10] ,[12]. Pero aho- 
ra, como nos propusimos, las variables de estado y los parámetros empíricos han 
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quedado determinados como funciones analíticas de los parámetros de nuestro 
modelo. 


(6)Discusión y conclusiones 


(a) Los resultados sobre la dinámica del umbral de excitación eléctrica de 
fibras nerviosas obtenidos en esta parte del trabajo permitieron relacionar 
las amplitudes umbral de corriente en los electrodos, en el caso de la esti- 
mulación eléctrica, y las pendientes umbral de las rampas de corriente en 
las bobinas, en el caso de la estimulación magnética, con parámetros elec- 
trofisiológicos de membrana básicos, con el radio de la fibra (a través de la 
constante de espacio homogeneizada), con el revestimiento de mielina y la 
extensión de los nodos de Ranvier, y con la distribución del campo eléctri- 
co paralelo a la fibra, producido por el sistema estimulador en el conductor 
de volumen de los tejidos (a través de las proyecciones del factor de forma 
y de la longitud del intervalo de influencia). Pese a que las fórmulas que 
dedujimos son aproximaciones analíticas, establecen relaciones entre va- 
riables que otros modelos no suministran, y en consecuencia podrían resul- 
tar de interés para el diseño de electrodos de implante o bobinas estimula- 
doras. 

Un problema para cuya solución cabe esperar que los resultados obtenidos 
en este trabajo, puedan aportar una perspectiva diferente a la que aportan 
otros enfoques, es el que se plantea a propósito de la estimulación diferen- 
cial de fibras nerviosas y musculares. 


(b)Las ecuaciones (35) y (36), junto con los grupos de fórmulas (33) y (34) 
obtenidas reformulando la linealización de la dinámica modal en el con- 
junto de decaimiento, junto con barrera umbral (barrera tal que una vez 
traspasada se origina un PA) no son otra cosa que una versión fundamen- 
tada de las teorías fenomenológicas clásicas de dos factores para la excita- 
ción mediante electrodos externos. Pero además obtuvimos una versión 
más general que la versión original, puesto que con las modificaciones del 
caso, resulta aplicable a la estimulación magnética. Las fórmulas (33) pre- 
dicen correctamente los resultados experimentales de los electrofisiólogos 
clásicos sobre las modificaciones que sufren los parámetros empíricos de 
la teoría de dos factores cuando se modifica la distancia entre el electrodo 
y la fibra, el tamaño y forma del electrodo, el radio de la fibra, y la concen- 
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tración externa de Ca** [4], [7], [8], [10] ,[12]. El estudio detallado de todos 
estos aspectos y de las implicaciones de las fórmulas analíticas obtenidas 
en este trabajo para las curvas intensidad duración, debería ser llevado a 
cabo en un tercer trabajo, al igual que los efectos debido a las variables de 
activación de la membrana. 


(c)Siempre que el electrodo o la bobina no se encuentren muy cerca de la 
fibra mielínica, cabe esperar que las fórmulas para las reobases, las cons- 
tantes de tiempo y los factores sigma (tanto para la excitación catódica de 
cierre y la excitación anódica de apertura), resulten aplicables para estudiar 
el efecto de la variación de los parámetros del sistema dispositivo estimu- 
lador-conductor de volumen-fibra blanco de la estimulación. La constante 


de tiempo depende de a a través de T yde A. Cuando — aumenta de0a 


1, la constante de tiempo varía desde la correspondiente a una envoltura 
completa de mielina hasta la de una fibra completamente desprovista de 
mielina. Con una adecuada selección de parámetros para la membrana, el 


caso P=1 se podría aplicar al estudio de la dinámica del umbral de una 


fibra muscular esquelética excitada directamente por el campo externo (es 
decir, no excitada a través de las sinapsis con fibras nerviosas). 

Los incrementos en la constante de tiempo y de la reobase con la longitud 
del intervalo de influencia que predicen las fórmulas obtenidas para la es- 
timulación eléctrica catódica y anódica, han podido ser confirmados en 
experimentos in-vivo, in-vitro e in-silico para el caso de estimulación eléc- 
trica con electrodos puntuales en condiciones en las que la longitud del 
intervalo de influencia es proporcional a la distancia entre el electrodo y la 
fibra [15] [16] [18]. Para estimulación magnética con equipos de tercera 
generación, al parecer todavía no hay datos adecuados. 


(d)La deducción de fórmulas para la curva intensidad-duración correspon- 
dientes a las excitaciones catódica de cierre y anódica de apertura se basó 
en la dinámica del umbral para las amplitudes del modo fundamental (co- 
rrespondientes al voltaje a través de la membrana y a la variable de recupe- 
ración), considerado como modo dominante y desacoplado de los demás 
modos. La excitación catódica de apertura queda fuera del alcance de esta 
simplificación. En la primera parte pudimos obtener y en parte estudiar la 
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excitación anódica de cierre añadiendo modos de orden superior al estudio 
de la dinámica del umbral. Algo análogo cabe esperar en relación con la 
excitación catódica de apertura, empleando el modelo que tiene en cuenta 
la variable de recuperación. 
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